Phase Transition in an Exactly Solvable Extinction Model 
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We introduce a model of biological evolution where species evolve in response to biotic interactions 
and a fluctuating environmental stress. The species may either become extinct or mutate to acquire 
a new fitness value when the effective stress level is greater than their individual fitness. The model 
exhibits a phase transition to a completely extinct phase as the environmental stress or the mutation 
rate is varied. We discuss the generic conditions for which this transition is continuous. The model is 
exactly solvable and the critical behavior is characterized by an unusual dynamic exponent z — 1/3. 
Apart from predicting large scale evolution, the model can be applied to understand the trends in 
the available fossil data. 



An important question that intrigues both scientists 
and non-scientists alike, is that of the origin and evolu- 
tion of life on this planet. Where did all the diverse vari- 
ety of species on earth come from? Why did some of the 
species which initially existed, got wiped out while oth- 
ers survived? What are the factors influencing the pro- 
cess of mass extinction? These have emerged as a multi- 
disciplinary field of research over the last few decades, 
attracting attention of researchers from various branches 
of science. The phenomenon of species extinction is an 
equally important event, and is an inextricable part of 
the evolution process. In recent times, several models of 
evolution and extinction have been proposed. Some early 
evolution models have considered evolutionary dynamics 
of interacting species on a rugged fitness landscape [11, Q ■ 
In these models, under repeated mutation and selection 
fit species tend to climb up on the fitness landscape un- 
til a local maximum or peak is reached. The landscape 
may also be coevolving with the evolution of the species. 
A pioneering work by Bak and Sneppen |3|] modeled the 
extinction of coevolving species on a fitness landscape as 
a self-organized critical (SOC) process. Here, coevolu- 
tion of species can trigger coevolutionary avalanches of 
extinction events and the avalanche size is distributed 
algebraically with an universal exponent r in the range 
1 < T < |[J]. However, a major drawback of these mod- 
els is that they completely ignore the role of the envi- 
ronmental stresses e.g., climatic, geological or exogenous 
stresses[l]. Recently, Newmanp has proposed a model 
considering only environmental stresses as the cause of 
large scale extinction, which can explain the trends in 
the available fossil records Q. 

In this article, we study a model of biological evolu- 
tion taking into account both external stresses as well as 
biotic interactions between species as contributory fac- 
tors for species extinction. Here, less fit species either 
become extinct or they mutate, due to changes in envi- 
ronmental stress level. Biotic interaction is incorporated 
as cooperativity among the species, which has not been 
considered in any of the models mentioned above. It 
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is well known that cooperativity is an important factor 
for proper functioning of every ecosystem which arises 
from the interdependence between species, e.g., via the 
food web^. Such interactions at the level of individual 
species are in fact important in modeling evolution and 
extinction on ecological time scales, as has been argued 
in Ref. @. 

We show that, a phase transition leading to a com- 
plete extinction of the species may occur in this model, 
as the environmental stress is increased. Generally such a 
transition occurs discontinuously and only under specific 
conditions it becomes continuous. For the continuous 
transition, the critical value of the stress, the the crit- 
ical exponents (a, /?, i^,,, z) and the scaling functions 
can be calculated analytically. The critical behavior of 
the system is found to be robust against the variation of 
mutation rate and the fluctuations in stress. 

Let us now describe the model in detail. We consider 
an ecosystem consisting of N different species and let 
Xi [i = l,2,...,iV) denote the fitness of the i^^ species. 
The fitness Xi is drawn from a fitness distribution J-{x). 
The species are subjected to a fluctuating environmen- 
tal stress S which is drawn from a distribution 'D{S). 
However, due to the presence of cooperativity, the ex- 
isting species in the ecosystem experience an effective 
stress s = S/C, where C is the measure of cooperativ- 
ity among the species. Thus, for a given environmental 
stress S, the effective stress is smaller when cooperativity 
among species is large and vice- versa. In general, C{Nt) 
depends on the existing number of species Nt at time 
t. If the number of existing species is very large, they 
compete for resources e.g., food, habitat etc. and their 
cooperativity decreases. Thus C{Nt) is expected to be a 
decreasing function for large Nt. However, for very small 
Nt, the species tend to cooperate for survival. In this 
regime, the resulting cooperativity C{Nt) is an increas- 
ing function. 

Extinction of species in our model occurs as follows. 
At any instant t, species whose fitness Xi is numerically 
smaller than the effective stress st = St/C{Nt), either be- 
come extinct with rate p, or mutate with rate (1 — p) to 
a new fitness Xi randomly drawn from the distribution 
F{x). Note that the renewed fitness value may be ei- 
ther same, greater, or smaller as compared to the current 
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value of the species fitness, and correspond to a neutral, 
favorable or an harmful mutation respectively. 

Since in our model the extinct species are not re- 
populated by new ones, the number of existing species 
Nt, and therefore, the density pt = Nt/N, can only 
decrease with time. The new effective stress is then 
St+i — Sf+i/C{Nt+i), where, 5*4+1 is a new environ- 
mental stress value drawn from The dynamics 
stops when either all the the surviving species are fit, 
or none of them survive under the applied stress. There- 
fore, depending on the environmental stress there is a 
possibility of complete extinction of the species. This 
transition, from a phase with finite population to com- 
plete extinction, may occur continuously or abruptly as 
(5) = / S'D{S)dS is increased. We will characterize this 
transition in detail considering a — (S) /N and p as con- 
trol parameters, and provide a phase diagram in the p - 
a plane. 

First, let us assume that, (i) the resources are infinite 
and therefore, the cooperativity C{Nt) increases mono- 
tonically, say linearly, with Nt, i.e., C{Nt) = Nt and (ii) 
the environmental stress S do not fiuctuate in time, i.e., 
T>{S) — 6{S — Sa). The more general cases, like different 
functional forms of C{Nt) and effects due to fiuctuation 
in the stress will be discussed at the end. 




FIG. 1. (Color online) Schematic description of the removal 
and renewal events. The a;-axis represents the cumulative 
distribution, which extends from g{Q) = to gipo) = 1. At 
a given time t, the species which are surviving (extinct) are 
denoted by solid (dotted) horizontal lines and a marker 774 = 
gicr / pt) is used to indicate the point below which the species 
are considered unfit. Of these surviving unfit species (their 
density is 74), 1—p fraction will mutate to have new fitness 
values extending over the whole range, increasing the density 
of surviving species at t -f 1 by (1 — p)7t. Effectively, at t + 1 
all the solid lines up to -qt are removed and one additional 
solid line with weight (1 — p)'yt is drawn. 



Starting from a initial density pQ = \, the density of 
existing species evolves as. 



where, 74 is the density of the unfit species at time t, 
i.e. N"ft is the number of existing species with fitness 
value smaller than Sa/Nt = (j/pt- Clearly, p fraction of 
A^7t species becomes extinct, and 1—p fraction undergo 
mutation, acquiring new fitness value spanning the entire 
range of distribution function J-(x). 

Effectively, the above dynamics amounts to the re- 
moval of all the Njt unfit species from the system at 
time t, and introducing 7V(1 — p)"ft species with renewed 
fitness value. The time sequence of such removal and re- 
newal events are schematically shown in Fig. [TJ The solid 
horizontal lines extending from 17(0) = to 5(00) = I rep- 
resent the density of species which underwent mutation 
at the previous time instant and is equal to (1 — p)-jt', 
g{x) = T{x')dx' is the cumulative distribution func- 
tion. The dashed lines correspond to the species which 
have been removed from the system. Utilizing this, it 
is easy to calculate the density of unfit species 7t+i and 
is equal to the weighted sum of the length of solid lines 
upto rjt, with the weight factors (1 -—p)"ft, where, 

/■o"/Pt 

?7t = / T{x)dx = gia/pt) 



as shown in Fig. ([T|). Thus, 

t-i 

i + {i-p)j2it' 



7t+i = ('7t+i - Vt) 



t'=0 



+ (1 - phtvt+i 

(2) 



Eqs. HI) and ([2]) describe the dynamical rules of our 
model and can be recast into a simple form using 



t 

E 

t'=0 



If, 



Pt+l Po 



TTt + l = ?7t[l + (1 -p)l^t\ 
pTTt+l = 1 -p77t[l + (I -_p)7rt] (3) 



1), this set of equa- 
?7(. This special case 



pt+i = Pt - Pit = Po - It' 

t'=0 



(1) 



When mutation is absent (p : 
tions Eq. ^ reduces to pt+i — 1 
of our model, with C{Nt)_ = Nt and V{S) = 6{S - Sa) 
has been studied earlier [!|] as the democratic Fiber Bun- 
dle Model (dFBM) in the context of failure processes. 
In the dFBM model, a heavy load weighing Sa hangs 
from a rigid anchor by a bundle of N elastic fibers, each 
having a certain breaking strength Xi. Initially, all the 
fibers are intact and share the applied load equally, each 
experiencing an effective load Sa/N. At each time step, 
weaker fibers (fibers with strength less than the effec- 
tive load) fail, and the load is re-shared equally among 
the remaining intact fibers. Thus the effective load per 
fiber increases creating an avalanche of failure events. If 
the initial load Sa is low, this process reaches a station- 
ary state with some intact fiber which eventually support 
the load, whereas a complete failure occurs for high Sa. 
Thus, at some critical value of Sa the dFBM model ex- 
hibits a breakdown transition which may be abrupt or 
continuous[10]. 

Coming back to the general case p ^ 1, let us first 
compute the fixed points of Eq. ([3]), where TTt+i = TTt ~ 
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TT* and pt+i — pt — p* ■ Thus, we have tt* = g{a/p*)[l + 
(1 — p)^*] and 

This equation may have multiple solutions for p*; the 
largest among them ps ~ Max(p*) is stable since the 
right hand side of Eq. Q is a non-decreasing function of 
p* bounded in the range (0, 1). Therefore, starting from 
the initial density po = 1, the density decreases and even- 
tually approaches a stationary value ps- 

The steady state density ps is the rightmost intersec- 
tion point of the curves 

l-5(l/x) 



y — ax and y = G(x) 



l-(l-p)g(l/x)' 



(5) 



where x replaces p*/(7 in Eq. ^ and G(x) is a non- 
decreasing function with G(0) = and G(oo) = 1. 
Clearly, these two curves intersect at x = for all values 
of a. If the curves intersect at other points (x > 0), 
the rightmost one, (say) x^ correspond to the steady 
state density ps = axg. This is described schematically 
in Fig. [5^. Existence of a solution Ps > indicates that 
the system is in a non-extinct phase while complete ex- 
tinction occurs when the only solution is ps = 0. 

The transition point between the two phases and the 
order of the transition can be determined from the behav- 
ior of G(x). Consider that G(x) is bounded from above 
by the a line y = ax such that G(x) < ax V x. Let this 
line be a tangent of G(x) at some x = x. If the line is 
a tangent at multiple points (see Fig. the rightmost 
one will be x. Clearly the system is in (non-) extinct phase 
for any a (smaller) greater than a. So, at the transition 
point ac = a the order parameter is ps = ax. The transi- 
tion will be discontinuous when x is nonzero, since in this 
case the order parameter has a nonzero value at the crit- 
ical point (see Fig.l^b). A continuous transition occurs if 
X = 0, and the critical point is, therefore, ac — a — G'(0). 
This implies that the transition is continuous only if the 
tangent line of G(x) at x = bounds the curve from 
above, i.e., when G(x)/x < G'(0) V x (see Fig.[2]a). 

It may be mentioned here that the extinction transition 
will be discontinuous if T{x) has a cut-off (say) at Xm, 
such that T{x > Xm) = 0. This implies that g{x > Xm) — 
1 and therefore, G(x < 1/xm) = 0. A typical form of 
such a function, with x^ = 1, is shown in Fig. [2l3(inset). 
In this case, the transition will occur discontinuously as 
i 7^ 0. 

In the following, we discuss the continuous phase tran- 
sition in detail. Since the critical point ac = G'(0), is 
related to f;(x) and its derivative as x — >■ cx), we expand 
g{x) as, 

g{x) = l~a{l/xf + ... (6) 

where a is a positive constant. Thus for A < 1, the critical 
point (Tc = G'(0) = oo and hence complete extinction 
can never happen. Again, for A > 1 the transition is 
discontinuous as in this case G'(0) — 0. Thus, the only 
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FIG. 2. (Color online) (a) A typical G(x), bounded from 
above by j/ = ax where a = G'(0), results in a continuous 
extinction transition with pa = ctxs. (b) For generic G(x) the 
transition is discontinuous; (inset) G(x) for a bounded fitness 
distribution (see text). 



non trivial extinction phase transition occurs for A = 1 
where the critical point ac — a/p is nonzero and finite. 
Note that the higher order terms in Eq. ^ do not play 
any role in determining the critical point. 

Let us now characterize the phase transition in terms 
of the critical exponents and the scaling functions. From 
Eq. (4), it is evident that for small e ~ ac — a, the order 
parameter ps = {e/ac)^ , with (5 = 1. The scaling func- 
tions can be derived from the dynamics of the model near 
the critical point, which can be approximated, following 
Eq. g]), by 

1 - .g(lM) 

= 1-(1-P).(1M); 
Taking the continuum limit and retaining terms upto 
second order in p, we have 

ialp~e)'^ = ep-'^^p\ (8) 
at 

where CTc = a/p. Close to the critical point (e 0), 
by rescaling the variable as g = e^^p and t = et, Eq. (|5]) 



can be written in a scaling form as -4r 

^ p at 



Q 



a (I-P) „2 

— 7? — Q ■ 



The formal solution of this rescaled differential equation 
can be written in one of the following forms; 



p{e,t) = 



(9) 



where, the exponents a = 1 = v^. From the scaling rela- 
tion P = ai^ii, we again get /3 = 1. The functions fa{x) 
and [x) can be obtained by solving Eq. ^ with bound- 
ary condition /9(e,0) = 1 which gives. 



fa{x) = Xffj{x) = 



9 

p X 



a'^{l — p) + exp{~'px/a) 



(10) 



A finite size scaling relation can also be derived for 
this model. For a system of size A'^, the mean number of 
existing species (Nt) will have a finite nonzero value even 
at the critical point ac (see Fig. [3]) , originating from the 
large critical fluctuations in Nf . Fluctuations of Nt result 
from the inherent stochasticity of the fitness values and it 
is expected that its width will be proportional to ^/{Ntj. 
These fluctuations can be incorporated in the model by 
modifying the density pt as. 
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(11) 



where, pt = {Nt) /N, is the mean density and fc is a 
constant of proportionahty. Now using Eq. in the 
Eq. ([7]) and keeping the leading terms in N and p, in the 

dp 



continuum hmit we have, — = k[y^ p/N ~ Kp^], which 

on a rescahng of variables g — N'^/'^p and t = tN^^^^ 
yields, 



(12) 



where K = a{l — p)/pk. The solution g{t,) can be 
rewritten in the scaling form 



p{N,t) = N-^/''^fi,{tN-') 



(13) 



where z = \ = (since (i 



1, we have v±_ = 3). 
Note that the dynamical exponent z satisfy the scaling 
relation z — v^.jv^. The unusually low dynamical expo- 
nent signifies the accelerated extinction occurring in this 
model. 
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FIG. 3. (Color online) (a) ps as a function of a for fitness 
distribution T{x) — exp(— l/a;)/a;^ and p — 0.8. Simulation 
results for different 7V(symbols) are compared with solution 
of Eq. Q (solid line), (inset) log- log plot of ps vs. N a.t ac ~ 
1/p is linear with slope — /J/i^x ~ —1/3. (b) Phase diagram of 
the model; complete extinction occurs for ct > CTc = 1/p. 

Until now, we have considered the system with infi- 
nite resources which however is unrealistic. For finite 
resources, the cooperativity function decreases with in- 
crease in the number of species Nt , and a generic form can 
be chosen to be C{Nt) = exp{-Nt/Nm), where the 
scale Nm is proportional to the total number of species 
N. In this case, Eq. ([IJ with p — 1 can be written 

as pt+i = 1 - g , where p,„ = N„i/N. 

\Pt exp(-pt/p,„)/ 



Clearly, in the thermodynamic limit, this dynamics will 
reach a stationary density < ps < 1 only for 6* = 1. 
The critical behavior in this case is identical to that of 
the system with infinite resources. 

To demonstrate the continuous extinction transition 
we choose a specific fitness distribution F{x) — x~^e"^/^ 
and a linear cooperativity function. The steady state 
density for p = 1 can be calculated exactly using Eq. 
dH) as Ps = 1 + (TlF(-cr~^ exp(-cr"^)) where W{x) is 
the Lambert W function. For p 7^ 1, the solution of 
Eq. Q is found to be in excellent agreement with ps ob- 
tained from numerical simulations (Fig. EJa)). The inset 
shows that p^ - N^f^/"^ with l3/vx_ = 1/3. The line of 
criticality a = 1/p separates the extinct phase from non- 
extinct one (Fig. Ellb)). This critical behavior is robust 
against fiuctuations which has been checked numerically 
by adding gaussian noise to the environmental stress. 

In conclusion, we have introduced an exactly solvable 
evolution model which incorporates some of the impor- 
tant features of biological evolution, namely, extinction of 
less fit species under environmental stress, cooperativity 
among the species and mutation. Under fairly general 
conditions, the model undergoes a discontinuous phase 
transition into a fully extinct state, whereas for certain 
specific choice of fitness distribution and cooperativity 
function the transition occurs continuously. The criti- 
cal point, the critical exponents (a, /3, i^,,, z) and the 
scaling functions of the continuous transition are calcu- 
lated analytically. In particular, the dynamical exponent 
z = ^ is quite unexpected and implies an accelerated 
species extinction. Apart from predicting large scale evo- 
lution via a critical dynamics [l3| and existence of a com- 
pletely extinct state, the model can also be applied to 
understand the trends in the available fossil data, as has 
been attempted by other models of evolution 4j . Close 
to criticality, the distribution of extinction events in this 
model follow a power-law with exponent t = | (obtained 
from simulation of the model). Since here the growth 
exponent ?^ = (no repopulation), one may argue [iT^ 



In fact, this value is consistent 



that T — -r: — — ■ — 

1+rj+a 2 

with the distribution of extinct families of marine ani- 
mal species, where r varies in the range 1.35 to l.QSflsf. 
Our model can be extended further to include other com- 
plex features e.g., repopulation of extinct species, local 
inter-species interactions in different dimensions, species 
dependent external stress, speciation. 

The authors thank M. Basu for providing help with the 
figures and U. Basu for comments on the manuscript. 



[1] S. Wright, Proc. Nat. Acad. Set. 58, 165 (1967). 
[2] S. A. Kauffman and S. Levin, J. Theor. Biol. 128, 1145 
(1987). 

[3] P. Bak and K. Sneppen, Phys. Rev. Lett. 71, 4083 (1993). 
[4] M. E. J. Newman and R. G. Palmer, adap-org/9908002 
[5] M. E. J. Newman, J. Theor. Biol. 189, 235 (1997). 



[6] J. J. Sepkoski, Milwaukee Public Museum Contributions 

in Biology and Geology 83 (1993). 
[7] B. Drossel and A. J. McKane, Handbook of Graphs and 

Networks (eds. S. Bornholdt and H. G. Schuster), pp. 

218-247. Wiley-VCH, Berlin. 
[8] D. Chowdhury and D. Stauffer, Phys. Rev. E 68, 041901 



5 



(2003). 

[9] S. Pradhan, A. Hansen, B. K. Chakrabarti, Rev. Mod. 

Phys. 82, 499 (2010). 
[10] J. V. Andersen, D. Sornette, and K. Leung, Phys. Rev. 

Lett. 78, 2140 (1997). 
[11] Solution of Eq. ((T2| is t = —^G^^Kq), where G{x) = 



2V3tan-(l±|^) + logit^. 
[12] M. A. Munoz, R. Dickman, A. Vespignani and S. Zapperi, 

Phys. Rev. E 59, 6175 (1999). 
[13] R. V. Sole and J. Bascompte, Proc. R. Soc. London B 

263, 161 (1996). 



